
#include "../../matrix/matrix.h"
#include "analisys.h"
/**
 * Compute the Global Otsu Threshold.
 * 
 * @param Hist Histogram
 * @return threshold value.
 */
real imp_image_get_otsu_threshold(ImpMatrix * Hist) {

	real * hist;
	real Th;
	long int sum;
	long int total;
	int n;

	hist = Hist->data;

	sum = 0;
	total = 0;
	for (n = 0; n < 256; n++) {
		sum += n * hist[n];
		total += hist[n];
	}

	long int wB, wF;
	long int sumB, sumF;
	float uB, uF;
	float var;
	float var_max;

	Th = 0;

	var_max = 0;
	wB = 0;
	sumB = 0;
	for (n = 0; n < 256; n++) {

		wB += hist[n];
		if (wB == 0) continue;

		wF = total - wB;
		if (wF == 0) break;

		sumB += n * hist[n];
		sumF = sum - sumB;

		uB = sumB / wB;
		uF = sumF / wF;

		var = wB * wF * (uB - uF)*(uB - uF);

		if (var > var_max) {
			var_max = var;
			Th = n;
		}
	}

	return Th;
}
real imp_image_get_otsu_threshold_from_image(ImpMatrix* Image) {

	real th;

	ImpMatrix * Hist = NULL;

	imp_image_histogram_mat(Image, &Hist);

	th = imp_image_get_otsu_threshold(Hist);

	imp_object_delete(Hist);
	
	return th;
}